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Abstract 

In this paper, we implement and test two types of market-based models for European-type options, based 
on the tangent Levy models proposed in p4i and O. As a result, we obtain a method for generating Monte 
Carlo samples of future paths of implied volatility surfaces. These paths and the surfaces themselves are free of 
arbitrage, and are constructed in a way that is consistent with the past and present values of implied volatility. We 
use a real market data to estimate the parameters of these models and conduct an empirical study, to compare the 
performance of market-based models with the performance of classical stochastic volatility models. We choose 
the problem of minimal-variance portfolio choice as a measure of model performance and compare the two 
tangent Levy models to SABR model. Our study demonstrates that the tangent Levy models do a much better job 
at finding a portfolio with smallest variance, their predictions for the variance are more reliable, and the portfolio 
weights are more stable. To the best of our knowledge, this is the first example of empirical analysis that provides 
a convincing evidence of the outperformance of the market-based models for European options using real market 
data. 


1 Introduction 

The existence of liquid markets for equity and volatility derivatives, as well as a well-developed over-the-counter 
market for exotic derivatives, generates a need for a modeling framework that is consistent across time and across 
financial instruments. Within this framework, once a model is chosen so that it matches both the present prices 
of liquid instruments and their past dynamics, it is expected to produce more realistic results for the problems of 
pricing and hedging of exotic instruments. In addition, such models can be used to quantify the risk embedded in 
portfolios of derivative contracts. Needless to say, evaluating and managing the risk of such portfolios is crucial for 
proper functioning of the financial markets: recall, for example, that VIX index, itself, is a portfolio of European 
options written on S&P 500. 

In this paper we investigate an arbitrage-free modeling framework for multiple European-type options written 
on the same underlying, which is consistent across time and products. In particular, this framework allows to 
resolve one of the nagging challenges of quant groups supporting equity trading: i.e. how to generate realistic 
Monte Carlo scenarios of implied volatility surfaces which are consistent with present and historical observations? 
As mentioned above, such models can be used to address the problems of pricing, hedging and risk management. 
Herein, we implement several such models using real market data and conduct a numerical experiment which 
demonstrates clearly the advantages of this modeling approach. 

The attempts to model the dynamics of implied volatility surface directly can be dated back as early as the 
“sticky smile model” and the “sticky delta model” (also known as “floating smile model”) (see Section 6.4 of 
1261 for the definitions). As an improvement of the two models, Cont et al. later proposed a multi-factor model 
of implied volatility surface in |[8l and ||9l, where they applied a Karhunen-Loeve decomposition on the daily 
variations of implied volatilities. It turns out that the first three eigenvectors could explain most of the daily 
variance, and a mean-reverting factor model based on the three eigenvectors is then constructed for future implied 
volatility surface. The major issue with these early attempts is that the proposed models for the dynamics of implied 
volatility are either too restrictive, not allowing to match the historical evolution of implied volatility, or too loose, 
so that they may contain arbitrage opportunities. While the importance of the first issue for any time-series analysis 
is very clear, the second one deserves a separate discussion. Indeed, what do we mean by arbitrage opportunities 
in a model for implied volatility and why do we need to avoid it? There are two types of arbitrage opportunities 
we refer to: static and dynamic. A given implied volatility surface contains static arbitrage if it is impossible to 
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obtain such a surface in any arbitrage-free model for the underlying. The fact that not every surface can be an 
arbitrage-free implied volatility simply follows from the well-known static no-arbitrage restrictions on the shape 
of a call price surface: e.g. monotonicity and convexity in strikes, etc (cf. im and ini). Notice that a violation 
of any of these conditions leads to an obvious arbitrage opportunity which is very easy to implement, hence, it 
is natural to assume that every implied volatility surface is free of static arbitrage. This, in turn, implies that any 
realistic simulation algorithm for future implied volatility surfaces has to produce surfaces that are arbitrage-free: 
otherwise, the algorithm generates outcomes that are simply impossible. The static no-arbitrage conditions are 
rather difficult to state explicitly, in terms of the implied volatility surface itself (without mapping it to a call or 
put price surface first). Nevertheless, it is not hard to deduce from the existing necessary conditions (cf. ED) that 
the set of arbitrage-free implied volatility surfaces forms a “thin” set in the space of all (regular enough) functions 
of two variables. Hence, it is a non-trivial task to construct a modeling framework that excludes static arbitrage 
in the implied volatility surface. The dynamic arbitrage adds to this problem, and it refers to a restriction on the 
evolution (i.e. the time increments) of implied volatility surface, rather than its values at a fixed moment in time. 
This restriction follows from the same arbitrage considerations for option prices. However, its associated arbitrage 
strategies are not as straightforward as in the case of static arbitrage. In addition, the simulated implied volatility 
surfaces that contain only dynamic arbitrage are, typically, very close to the ones that are arbitrage-free, when the 
time horizon is small (it is related to the fact that dynamic arbitrage only changes the drift term of the implied 
volatility, which is much smaller than the diffusion term, for small times). This is why, eliminating the dynamic 
arbitrage in a model for implied volatility surface is often viewed as a “second priority” for risk management. 
Nevertheless, we believe that a good model should exclude both types of arbitrage, in order to produce realistic 
dynamics of implied volatility surface (for risk management) and eliminate the possible arbitrage opportunities 
(for pricing). 

We have already mentioned that it is not a trivial task to construct a model of implied volatility that excludes 
arbitrage opportunities. In fact, when trying to model the surface directly, the first challenge that one faces is: 
how to describe the space possible implied volatility surfaces? Note that, as discussed above, the existing char¬ 
acterizations of arbitrage-free implied volatility surfaces are rather implicit. In addition, if the resulting space is 
not an open subset of any linear space (which it is not), what kind of mathematical tools can be used to describe 
evolution in space? Recall, for example, that all statistical models of time-series are defined on linear spaces (or 
those that can be easily mapped in to a linear space). Hence, it appears natural to map the space of possible implied 
volatility surfaces to an open set in a linear space, and then proceed with the construction of arbitrage free models. 
Such mapping became known as a code-book mapping, and it turns out that it can be constructed by means of 
the so-called tangent models (cf. El, a, 0). The concept of a tangent model is very close to the method of 
calibrating a model for underlying to the target derivatives’ prices prices (in the present case, European options 
calls). Consider a family of arbitrage-free models for the underlying, M.{9), parameterized by 6, taking values 
in a “convenient” set 0 (an open set of a linear space). For any given surface of option prices (or, equivalently, 
any given implied volatility surface), we can try to calibrate a model for this family to a given surface of option 
prices (or, equivalently, to a given implied volatility surface). In, other words, we attempt to find 6 G Q such 
that: C^{T, K) = C{T, K), for all given maturities T and strikes K, where C{T, K) is the given call price, and 
C^{T, K) is the call price produced by the model M.{9). If the above calibration problem has a unique solution, 
we obtain a one-to-one correspondence between the call price surfaces and the models in a chosen family: 9 -n- C^. 
For every call price surface C = , the associated (calibrated) model Ai{9) is called a tangent modelj^ Notice 

that is always arbitrage-free, hence, we obtain the desired code-book mapping C = i-G 9. Now, the problem 
of static arbitrage has been resolved, and one simply needs to prescribe the distribution of a stochastic process {9t), 
taking values in a convenient set 0, in order to obtain a model for the dynamics of call prices {Ct = C^*), and, 
in turn, the dynamics of implied volatility surface. Finally, one needs to characterize all possible dynamics of {9t) 
that produce no dynamic arbitrage in the associated call prices (C'®*)- An interested reader is referred to 0, for 
a more detailed description of this general algorithm, and, for example, to ||2|, ffl, iflTl . 1^ . Il22l . Il27ll . for the 
analysis of specific choices of the families of models {M.{9)}. 

The idea of modeling prices of derivative contracts directly dates back to the work of Heath, Jarrow and Morton 
M, who analyzed the dynamic of bond prices along with the short interest rate. Such models have become known 
as the market-based models (or simply market models), as opposed to the classical spot models, since the former 
are designed to capture the evolution of the entire market, including the liquid derivatives. This approach has been 
extended to more general mathematical settings, as well as to other derivatives’ markets. The list of relevant works 
includes im, ED, ES, ED, Ea, lEl, in addition to those mentioned in the previous paragraph. Even though 

* It is important to remember that any such model serves only as a static description of option prices, and it does not describe their dynamics! 
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the notions of code-book and tangent models never appear in these papers, almost all of them follow the algorithm 
outlined in the previous paragraph (and described in more detail in ||3l), in order to construct a market-based model. 

Even though various code-books for implied volatility surface (or, equivalently, for call price surface) have been 
proposed and the corresponding arbitrage-free dynamics have been characterized, it was not until very recently that 
some of these models were implemented numerically. As is shown in the rest of the paper, the lack of such results 
is not a surprise given the complexity of the models. So far, the numerical implementations are mostly based 
on tangent Levy models proposed in il and 0: as the name suggests, this corresponds to a code-book which is 
constructed using non-homogeneous Levy (or, additive) models as the tangent models. Karlsson M implements a 
class of tangent Levy models with absolutely continuous Levy densities and no continuous martingale component. 
Zhao ll^ and Leclercq f22\ . on contrary, implemented the tangent Levy models whose Levy measure is purely 
atomic in the space variable. As opposed to ||33l, the work of Leclercq allows for tangent models with 
continuous martingale component and includes options with multiple maturities, but it does require that the Levy 
density possess certain symmetry, which may limit the ability of the model to capture the skew of the implied 
smile. All of the works lITSl . ll^ . estimate the parameters of the model from real market data. In addition, 
II 22 I conducts a numerical experiment comparing the performance of a market-based model to a classical spot 
model. The actual results of this experiment, however, do not provide a convincing evidence in favor of the 
market-based approach. We believe that the latter is simply due the choice of experiment and to the dehciency of 
the theory, and we intend to demonstrate it in the present work. 

The purpose of this paper is to propose implementation methods for two classes of tangent Levy models - 
with continuous and discrete Levy measures. These methods provide practical algorithms for simulating future 
arbitrage-free implied volatility surfaces, which are consistent with both present and past observations. Our hrst 
method is similar to the one used in ifTSl . but with a different “dynamic htting” part, and the second method is 
in the spirit of ll22ll . although we avoid the assumption of symmetry of a Levy measure made in ll22l . However, 
the most important original contribution of this paper is the numerical experiment which uses real market data to 
demonstrate clearly the advantages of market-based models for implied volatility (or, option prices), as compared 
to the classical spot models. To the best of our knowledge, this is the hrst convincing empirical analysis that 
justihes the use of market-based approach for modeling implied volatility surface. 

The rest of the paper is organized as follows. Section 2 starts by reviewing the work on tangent Levy models 
with continuous Levy density and continuous martingale component, developed in 13]. Then, we introduce the 
implementation approach for this models, which is based on double exponential jump processes, hence the name 
“Double Exponential Tangent Levy Models”. Section 3 introduces the implementation method for tangent Levy 
models with discrete Levy density, called “Discrete Tangent Levy Models”. The two approaches are then tested 
against a popular classical model in a portfolio optimization problem in Section 4. Section 5 concludes the paper by 
highlighting the main contributions and the future work. Appendices A-C contain technical proofs and derivations. 
Appendix D contains all tables and graphs. 


2 Double exponential tangent Levy models 

2.1 Model setup and consistency conditions 

In this subsection, we review and update the results of a, which serve as a foundation for the analysis in sub¬ 
sequent sections. Herein, we assume that the interest and dividend rates for the underlying asset are zero. In the 
implementation that follows, we discount the market data accordingly, to comply with this assumption. As in a, 
we denote by {St)t>o a stochastic process representing the underlying price, and assume that the true dynamics of 
S under the pricing measure Q are given by: 

St = Sof f Su-{e^ — i)[M{dx,du) — Ku{x)dxdu]. (2.1) 

Jo Jr 

Here, M is a general integer-valued random measure (not necessarily a Poisson measure!), whose compensator is 
Ku,uj{x)dxdu, where {Kt)t>o is a predictable stochastic process taking values in the function space Bq, dehned 
in ( |6.1| l. 

Lor any fixed time t > 0 and a given value of St, a stochastic process (iST)T>t is said to be tangent to the 
true model {St)t>o if the time-f prices of all European call options written on S can be obtained by pretending the 
future risk-neutral evolution of the index value is instead given by {ST)T>t from t on. Throughout this section. 
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for any fixed f > 0, we assume that the tangent processes S is in the form 

St = St + / Su-{e^ — [Nt{dx, du) — Kt{u, x)dxdu], 


( 2 . 2 ) 


Jt Jb 

for T G [0, T], where Nt{dx, du) is a Poisson random measure associated with the jumps of log S whose com¬ 
pensator is given by a deterministic measure Kt{u,x)dxdu. Notice that the law of S is uniquely determined by 
{St, Kt). Let Ct*’'^^ {T, x) denote the option prices generated by {Su)u>t, i.e 


:=E (5t - e")+|5t = 


VT > t, X G 


(2.3) 


The concept of a tangent model, then, requires that, for each fixed t G [0, T), 

Cf^’'‘^(T,x) = E [{St - e^)+\J^t] , VT > f, Vx e M. (2.4) 

Thus, at each time t, we obtain the code-book for call prices, given by {St, Kt)- Of course, the value of the code¬ 
book may be different at a different time t. Hence, we consider the dynamic tangent Levy models characterized 
by a pair of stochastic processes {St, Kt)tG[o,T] that satisfies ( |2.4| ). Here, S' is a positive martingale with dynamics 
given by (HB; K is progressively measurable positive stochastic process taking values in B (cf. ( |6.2| l). The 
dynamics of St and Kt are given by 


St = So-G /g Su- (e'^ - 1) [M{dx, du) - Ku {x)dxdu], 
Kt{T, x) = ko{T, x) + f* au{T, x)du + Jo K{T, x)dB^, 


(2.5) 


where (atjjgjp j.] is a progressively measurable integrable stochastic process with values in B, and, for each n G 
{1, • • • , to}, (/3[‘)tg[o T] ^ progressively measurable square integrable stochastic process taking values in T-L (cf. 
1 ^). 

Notice that p.5| l defines the dynamics of the code-book {St, Kt)t£^[o,TG hut it does not ensure that it does, 
indeed, produce tangent models at each time t: in other words, there is no guarantee that ( |2.4| i holds. Thus, 
additional “consistency” conditions have to be enforced to obtain models which are, indeed, tangent to the true 
underlying process. As shown in a, this consistency is, in fact, equivalent to the fact that call prices generated by 
these tangent models are free of dynamic arbitrage. In order to present the main consistency result, we state the 
following regularity assumptions on /3. 

Assumption 1. For each n < m, almost surely, for almost every t G [0, T], we have: 

RAl s}xpTe[t,T] ill l/3"{T,x)ldx < oo 

RA2 For every T G [t, T], the function (5f{T, ■) is absolutely continuous on K \ {0}. 

RA3 For any T G [t, f], (e“ - 1) P]^{T, x) = 0. 


Finally, we introduce some extra notation and formulate the consistency result, which is a simple corollary of 
Theorem 12 in 0. 


PtiT,x) := / j3'f{u,x)du. 


( 2 . 6 ) 


/tAT 


Theorem 1. ( Carmona-Nadtochiy 2012) Assume that {St)t^\^o T] martingale, jd satisfies the above regu¬ 

larity assumptions RA1-RA4, and Kt{T, x) > 0, almost surely for all t G [0, T) and almost all {T, x) G [t, T] x M. 
Then the processes {St, Kt)t^[o T] ^otisfying (2.5i are consistent, in the sense that ( |2.4| ) holds, if and only if the 
following conditions hold almost surely for almost every x S M and t G [0, T), and all T G {t, T]: 


1. Drift restriction: 


at{T,x) 


E 


Pl{T,y)pnT,x-y)dy 


PnT,x)- / P^{T,z)dz-P^{T,x)- / ~p'^{T,z)dz\. 


(2.7) 


4 





2. Compensator specification: Kfix) = Kt{t,x). 


Theorem[^ along with equations ( |2.5| l provide a general method for constructing a market-based model for call 
prices (i.e. an arbitrage-free dyanimc model for implied volatility surface). Indeed, choosing (f]} ..., /^™)tg[o t], 
we use the drift restriction in Theorem and the second equation in ( |2.5| l to generate the paths of (Ki)tG[o t]- 
Finally, to generate the paths of (S't)jg[o t]’ use the compensator specification in Theorem|^and the hrst 

equation in ( |2.5| ), after representing the random measure M through its compensator K and a Poisson random 
measure N (as shown in Q). However, in the present paper we avoid simulating (S't)tG[oT] all, by simply 
noticing that 




St .fit / 


{T,x + \ogSt)=E {ST/St-e^)+\St = St 


= E 


{ST-en^\St = l 


X + log St; a) = x; a), 

where x) is the Black-Scholes price at time f of a call option with maturity T and strike given that the 

level of underlying is at St and the volatility is cr. At any time t, regardless of the value of St, if we find the level 
of a that makes the right hand sides of the two equations above coincide, then the option prices in the left hand 
sides have to coincide as well. This means that we can obtain the implied volatility of in the maturity and 

log-moneyness variables, by computing the corresponding implied volatility of for which we do not need to 

generate St- 


2.2 Implied volatility simulation with tangent Levy models 

We first introduce the general framework of the simulation procedure. Our procedure has two stages, estimation 
and simulation. The estimation stage, where the additive density of the tangent process as well as its dynamics are 
fitted to market data, is performed in two steps: 

• Static fitting. In static fitting, the additive density Kt for each day t is obtained by least squares optimization 
which minimizes the squared difference between model prices and actual market prices. Notice that for any 
given day t, Kt is fixed and there is no dynamics involved, which explains the term ‘static’. 

• Dynamic fitting. In dynamic fitting, we recover the dynamics of the time series (kj). In view of the drift 

restriction in Theorem this boils down to determining the volatility terms This is done by 

applying the Principle Components Analysis to the time series of {Kt)t- 

Once the estimation is completed, we generate the future paths of {Kt) using Euler scheme Monte Carlo applied 
to the second equation in ( |2.5| l. From the simulated additive densities, we compute call prices and, then, 
implied volatilities by inverting the Black-Scholes formula. 

Within the general framework, the simulation stage is generic, but the static part of the estimation stage can 
be quite different depending on the specific subclass of tangent Levy densities k{u, x) that we ht to option price 
at any given time. In this section, we implement the procedure with the Levy densities arising from the double 
exponential Levy models proposed by Kou in 1191 . The small number of parameters in double exponential models 
and the availability of an analytical pricing formula for call options make the resulting family of tangent Levy 
models fairly easy to calibrate. 


2.3 Market data 

We use SPX (S&P 500) call option prices provided by OptionMetrics, an option database containing historical 
prices of options and their underlying instruments. Throughout the paper, we use the option data from two time 
periods: Jan. 2007 - Aug. 2008 and Jan. 2011 - Dec. 2012. Table[T]gives a quick summary of the two periods. We 
cut off the first period at August 2008 to reduce the impact of the financial crisis. 

On each day of a period, we only keep the options with time to maturity less than one year, whose best closing 
bid price and best closing offer price are both available, and take the average of the two prices as the option price. 
To ensure the validity of all prices, the contracts with zero open interest are excluded. As a result, there are roughly 
10 to 80 call contracts with valid prices available for each maturity. The log-moneyness (more precisely, the put 
log-moneyess, dehned as log(iT/S't)) of these call options ranges roughly from -0.3 to 0.1, varying for different t 
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and T. Our calibration also requires dividend and interest rate data available on OptionMetrics and the homepage 
of U.S. Department of Treasury, respectively. This dividend yield is recovered from option prices via put-call parity 
with the method proposed in Hi. On day t, we denote the dividend yield by qt, and the risk-free rate between t and 
T by Tt T- To simplify our implementation, we perform a simple transformation on the market data so that we can 
assume that the interest and dividend rates are both zero from now on; 

= with x = x - - qt)iT - t), (2.8) 

where C'^^'^{T,x) is the market price of a call option with maturity T and strike e“. The adjusted call prices 
CY^^*{T,x), corresponding to maturity T and strike e®, are then consistent with the assumption of zero interest 
and dividend rates (i.e. they do not contain arbitrage under thee assumptions). IN a similar way, we define the 
adjusted bid and ask prices, and . 

In this Section and Section]^ we will perform the calibration of the tangent Levy models on the time span from 
Jan. 3, 2007 to Dec. 31, 2007, denoted by [fg, T], In Sectionj^ data from both periods will be used to test the 
performance of the tangent Levy models. 

2.4 Static fitting 

Before we proceed with the static fitting, let us first have a quick review of the double exponential model. In such a 
model, the logarithm of underlying follows a pure jump Levy process whose jump sizes have a double exponential 
distribution. More specifically, assuming no diffusion term, the dynamics of the underlying are given by 


dSt = IJ. St-dt -\- St- d 


Nt 


^(exp(y,) - 1) 


(2.9) 


where p, is the drift term, Nt is a Poisson process with rate A, {Yt} is a sequence of i.i.d. random variables with 
asymmetric double exponential distribution, independent of Nt- The density of an asymmetric double exponential 
distribution is given by 

fviy) = P ■ Xie-^^yiy>o + q ■ X 2 e^^yiy<o, (2.10) 

where p,q > 0, p + q = 1 represent the probabilities of positive and negative jumps, and Ai > 1, A 2 > 0 are the 
parameters of the two exponential distributions. In other words, a double exponential model is a martingale model 
for the underlying whose logarithm is a pure jump Levy process, with the Levy density 

r]{x) = X{p ■ Aie"^i^l3;>o + q ■ A 2 e^"^la;<o). (2.11) 


One of the advantages of double exponential models is the availability of analytical pricing formulas for Euro¬ 
pean call options, which could greatly simplify the calibration. ifT^ gives the pricing formula for double exponen¬ 
tial models with a diffusion term. A minor modification of the derivation in ifT^ gives us the pricing formula in 
absence of the diffusion term, as shown in the lemma below (its proof is given in Appendix B). 

Lemma 1. Under the assumptions of zero interest and dividend rates, assume, in addition, that the underlying 
process S follows a double exponential process with Levy density given by ( |2.11[ ), under the risk-neutral probability 
measure. Then, the price of a European call option with strike K and maturity T is given by 

logiT) = St^ (-AC, A*,P*, At, a;; log {K/St) ,T-t) 

- (-AC, X,p, Ai, A 2 ; log (K/St) ,T-t), (2.12) 


where 


P Ai 

1 + c ■ Ai-r 


A* = A(C+1), C 


At — Ai — 1, At — A 2 + 1, 

pXi qX2 _ 

■ Ai - 1 A 2 + 1 
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and the function 'P is given by: 


with 


and 


A2;a,T) 


oo n 


= TTola-AiT^O + ^ ^ Pn,k 


n—1 k—1 


d'T)y Xi(a-^T)-, I 

-7]-e aa-^T>0 + 


la—/xT <0 


. 2=0 


£ i: f 1 - E !___ 


fiT<0^ 


n—1 k—1 


2=0 


Pn,k = Y. 


i—k 

n—1 


n — k — 1\ fn 


i — k 


Qn,k — ^ ^ 
i—k 

Pn,n = P^, Qn.n = ■ 


i — k 
n — k — 1\ fn 


e-^'^{XTy 

n\ 


s. \ i—k / \ \ n—i 

M \ { M 


Ai + A 2 

Ai 

Ai + A 2 


Ai + A 2 


n-i / , \ i-k 

M 


Al + Ag 


p^q^-\ l<k<n-l, 


l<k<n-l, 


(2.13) 


For each T/, with I = 1,..., L, we would like to find the set of parameters {A, Ai, A 2 ,p} that minimizes the 
difference between the market and the model prices. For practical reasons, we will work with time values instead 
of options prices. The market time value and the model time value are calculated as follows 




mkt^j 


K 


>^XlX2,P,j 


(TO = Cr"‘(r,,e"0-(^t-e"£+, 

(7)) = - [St - e"£+. 


There are two reasons for working with time values. Firstly, the time values go to zero for very large and very 
small log-moneyness, which allows us to truncate the x-space with negligible numerical errors. Secondly, time 
values and option prices are often of different magnitudes, especially for in the money options, with option prices 
much greater than time values, hence, working with time values is likely to result in smaller numerical errors. For 
fixed time t and fixed maturity Ti, the optimization problem can be written as 


min 

A>0,Ai>0.A2>0,pG(0.1) 


N 

i=i 


^j\Vt 


A,At .Ao 




(2.14) 


where ujj = \ (Ti, (Ti, ) | ' are the weights we put on different options to take into account 

the difference in liquidity (measured by bid-ask spread). For every fixed maturity Ti, the solution of the above 
optimization problem, (A*, A£ A£p*), yields the Levy density r]t{Ti,x) via ( |2.1 l| i. Then, we search for a function 
Kt{-, ■), such that 

1 fP 

r]t{Ti,x) = J Kt{u,x)du, (2.15) 


for every maturity Ti and all x S K. The resulting tangent model on day t is defined as a martingale model for 
the underlying whose logarithm is a pure jump additive (non-homogeneous Levy) process, with the Levy density 
Kt{-, •). It is easy to see that the call prices produced by this model, for every maturity Ti and strike coincide 


A^ A^ A^ 

with the prices produced by the double exponential model, Cj ’ {Ti, e^^ ). Thus, for a given t, the problem 

of static fitting is essentially a series of optimization problems ( |2.14| i, over all maturities Ti, along with the fitting 
problem ( |2.15| l. 

At the first glance, the optimization in ( |2.14| ) seems to have four parameters. However, the following constraints 
will reduce the number of parameters to two in our calibration; 
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• To improve the stability of small-jumps intensity over time, we would like the Levy density p(T;, x) to be 
continuous in x. The continuity at a; = 0 requires 


p • Ai = (1 - p) • A2 A2 


p 

1 -p 


Ai. 


(2.16) 


• In view of the results in Section 2.1 we have to impose the symmetry condition RA3 on /3”’s. A simple 
application of Ito’s lemma shows that, for the symmetry condition RA3 to hold, it suffices to choose every 
Kt, SO that 


J (e^ — l)Kt(T,x)dx 


is a deterministic function of T — f, for all times 0 < t < T <T. To achieve this, in view of ( |2.15| ), we need 
to choose every r]t{Ti, ■) so that the symmetry index 


- (T -t) ■= [ - l)r]t{T, x)dx = A ^ ^ ^ 

Jr \ Ai — i A2 + i 

is a deterministic function ofT — t. This yields; 

-(l + S(r-f)/A)(Ai-l) 

^ S(r-f)/A(Ai-l)2-2(Ai-l)-l’ 

where S is a hxed (estimated a priori) function. 


(2.17) 


(2.18) 


With the two constraints, our calibration takes only two variables: A and Ai. The condition p € (0,1) trans¬ 
forms to the following condition on Ai: 


r ( 1 , 00 ), if E{T-t)<0, 

+ if S(r-f)>0. 

As a result, the optimization problem ( |2.14| l can be re-written as 

A>U,AiEiAi . 


(2.19) 


( 2 . 20 ) 


where is the interval defined in \2.\9\ . The symmetry index function S(t), for all r G IR_|_, can be obtained on 
the hrst calibration day t = 0, solving a three-variable optimization problem. 


min 

A>0,Ai>l.pG(0,l) 


N 




( 2 . 21 ) 


and setting 


m) = A 


p 

Ai-l 


1 -P 
A 2 -f 1 


( 2 . 22 ) 


for every maturity Ti, and, finally, interpolating linearly between every T/_i and Ti. We summarize the calibration 
procedure for {pt(T;, •)} in the following algorithm: 


Algorithm 1: Algorithm for calibrating {T]t{Ti,-)} 
r Preprocess the market data according to ( |2.8| l; 

2 For f = 0, tTin the three-variable optimization p.21| ), without the symmetry condition, for all maturities, and 
compute S(.) by ( |2.22| i and linear interpolation; 

3 For the subsequent days t G (0, f], run the two-variable optimization p.20| i, with already estimated S, to 
obtain the time series of Levy densities (pt)te[o fj- 


Below are the calibration results. The Levy densities p on Jan. 3, 2007 - the first day of calibration - is obtained 
by the three-variable optimization (|2.21|l. From the calibrated parameters, we compute the symmetry index S via 
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( |2.22[ ), which is shown in Figure With the symmetry index S, we run the two-variable optimization ( |2.20| l on 
the following day, Jan. 4, 2007, and obtain the Levy densities p shown in Figure]^ The corresponding time values 
are shown in Figure We can see that the calibration results are quite precise in the sense that the time value 
falls between the bid and the ask values most of the time. As for the calibrated Levy densities ry, its values tend 
to decrease as the time to maturity increases (cf. Figure]^. The magnitude of S ( which measures the “asymetry” 
of the Levy measure) is decreasing with maturity as well. Both results are in line with empirical hndings on jump 
intensities and volatility skews. 

Next, for every day t, we need to hnd Kt that satishes ( |2.151 l. Notice that, if rit(T, x) is differentiable in T, we 
obtain; 

Pt(T,x) + (T-f)^^^^p^ = K*(r,x), (2.23) 

for each a; G M. The relationship ( |2.23| ) can be used to back out the additive densities jq t] from the calibrated 
Levy densities {r]t)te[o,T]- However, the calibrated densities r]t{T, •) are only dehned for T = Ti, hence, we need 
to interpolate them across maturities. An analysis of the calibrated Levy densities shows that rit{T,x) generally 
exhibits one of the following two patterns as a function of T. 

• For small jump sizes x, r]t{T, x) decreases rapidly as T increases. To ensure that the recovered k is non¬ 


negative, we used a combination of exponential function and power function 

r]t{T,x) = Ci{T - -f C 3 (r - t) exp(-C 4 (T - t)) -f cs (2.24) 

to ht 77, for any hxed x. The corresponding Levy density k can then be computed as 

Kt{T,x) = Ci(c 2 + 1)(T - + exp(-C 4 (T - t))(2c3(T - t) - csc^T - f)^) -f C5. (2.25) 

• For large jump sizes x, r]t{T, x) increases as T increases. The function we used to fit this scenario is a simple 
polynomial function 

rit{T, x) = Ci(r - ty + C 2 (r - tf + C 3 (r - tf + c^T -t)+ Cg. (2.26) 

Then, k is computed as 

Kt(r, x) = 5ci(r - ty + 4c2(T - ty + 3c3(T - ty + 2c4(T -t) + eg. (2.27) 


An illustration of the two scenarios together with an example of the reconstructed k is shown in Figure]^ 

2.5 Dynamic fitting 

Recall that, in view of ( |2.5| l, the Levy density n has the following dynamics: 

m 

dKt{T, x) = at{T, x)dt + ^ x)dB'^. (2.28) 

n—1 

In the dynamic htting, we need to assume that the time increments of k are stationary, which is only natural if we 
work with the time to maturity t = T — t instead of the maturity T. Namely, we dehne kt{T, x) = Kt{t + t, x) 
and its dynamics 

m 

dkt{T, x) = at{T, x)dt + P'tij, x)dBy (2.29) 

n—1 

A simple application of Ito’s formula shows that 

at(r, x) = at{t + r, x) + ^nd /3”(r, a;) =/3"(f-f r, a:). (2.30) 

To simulate future implied volatility surfaces, all we need are the diffusion terms /3"’s, because the drift term a 
can be computed from ^"’s. We assume that /3 ”(t, x)’s are deterministic and constant as functions of t, for any 
(r, x) (from a hnite family of points). Then, every increment Akt = kt — kt-i is a sum of a Gaussian random 
vector, corresponding to the diffusion part, and a vector that corresponds to the drift term (we view every surface 
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as a vector whose entries correspond to different values of (r, x)). Notice that the distribution of the Gaussian 
component is completely determined by its covariance matrix, hence, we will aim to choose /3"’s to match the 
estimated covariance matrix. Assuming that the drift term is bounded, it is easy to notice that the standard estimate 
of the covariance of Akt also provides a consistent estimate of the covariance of the aforementioned Gaussian 
vector, asymptotically, as the length of the time increments converges to zero. In the actual computations, we use 
daily increments - these are small compared to the time span of the entire sample, which is one year. To fit /3"’s 
to the estimated covariance matrix, it is natural to use the Principal Component Analysis (PCA), which finds the 
directions that explain most of the variance in the increments Akt- However, the PCA can not be applied directly 
because the number of points on the surface is close to the sample size, which is 251: for each t, we have call 
prices for 10 maturities and 21 jump sizes, which gives us 210 points on the k surface after static fitting. To reduce 
the number of points, we pick every other maturity and the 7 jump sizes whose intensities are larger than others 
across time t. This gives us 5 * 7 = 35 points on the reduced surface of {AK(}jg[g 

Applying PCA on the reduced surface, we see that the first three eigenmodes {/"(t, explain over 93% 

of the daily variance of k, as shown in Figure|^ To extend the values of the eigenmodes to other points (i.e. other 
jump sizes and maturities), we simply perform a linear interpolation. The first three eigenmodes have very unique 
characteristics. The first eigenmode takes the most prominent feature of k - the densities are concentrated around 
small jumps at very short time to maturity. This eigenmode can be understood as a combination of the “level” 
factor and the “slope” factor (appearing in a typical PCA result for yield curve dynamics) along both the maturity 
and the jump size directions. The second eigenmode shows the curvature along the jump size direction, and the 
third eigenmode shows the curvature along the time to maturity direction. As the eigenmodes {/"(r, are 

normalized, to obtain the diffusion terms /3"’s, we need to multiple the eigenmodes by the loading factors: 

^tiT,x) = ^/Xn-P{T,x), n = l,2,3. (2.31) 

Once we have /3"’s, we change the variables to pass to /3”’s and calculate the drift term a according to 
Figurej^shows the drift term a computed according to ( |2.7| i. Notice that a can then be computed as 

- / \ t , dKt{t + T,x) 

at(T,x} = at{t + T,x) -\ - — -, (2.32) 

where we have no problem with evaluating the partial derivative, as, in the static fitting stage, Kt was interpolated 
across maturities. 


2.6 Monte Carlo simulation of implied volatility surfaces 

Once all the terms in the right hand side of ( |2.29| ) are estimated, we can, for example, apply and explicit Eu¬ 
ler scheme to simulate the future Levy densities kt. However, we need to ensure that the simulated k^’s stay 
nonnegative at all times. Inspired by a, we incorporate a scaling factor in ( |2.29[ ) as follows: 

m 

dkt{T,x) = ^t^t{T,x)dt + 7t X! (2.33) 

n—1 


where 



inf kt{T,x)Ae 

tG [0,t] 


(2.34) 


with e = le“® and f = 1. Of course, this modification changes the diffusion term of kt, which was estimated from 
historical data. However, the value of e is chosen to be so small that, in the historical sample, 74 is always equal to 
one. Hence, if we use the /3"’s chosen in the previous subsection, the resulting dynamics are still consistent with 
the past observations. It is also easy to see that, since 74 is a scalar, the drift restriction ( |2.7| ) is satisfied by the new 
drift and volatility of k. Finally, this modification ensures that kt is almost surely nonnegative for any t. 

To simulate future values of k, we apply the explicit Euler scheme to (|2.331l, to obtain 


kt+At{T,x) = kt{T,x) +jfatiT,x)At + 'yt'^P''{T,x)ABt, (2.35) 

n—1 


with At being one day. Having simulated kt, we compute r/t via ( |2.15| ). Then, for every fixed maturity T, 
the option prices in the model given by the Levy density 774 (T, •) can then be computed, for example, using the 
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methods proposed in m or M- These methods are based on Fourier transform and can be implemented efficiently 
via numerical integration]^ In particular, in our simulation, we use the following formula to calculate future option 
prices: 






du 


Re 


exp {—iux) ipt \T,u 



(2.36) 


where (pt is the characteristic function of an exponential Levy process with the Levy density r]t{T, •), starting from 
one: 


(t)t{T,u) = exp 


-iu{T-t) / r)t{T,x){e^-l)dx + {T-t) / - l)r)t{T, x)dx 


From the above option price, {T,x), we can easily calculate the implied volatility by inverting the Black 
Scholes formula, assuming that = 1 and the interest and dividend rates are zero. As discussed at the very end of 
Subsection |2.1| this value is the same as the value of implied volatility of a call option for spot level St, strike Ste^, 
and maturity T, regardless of what the level of St is (hence, we don’t need to simulate it). Using this method, we 
simulate the implied volatility surfaces five days into the future starting from Dec. 13, 2007, as shown in Figures]^ 
and ID 


3 Discrete tangent Levy models 

3.1 Model setup and consistency conditions 

In this section, we work with a different class of tangent Levy models in order to solve the same problem - develop 
a consistent Monte Carlo simulation algorithm for the future implied volatility surfaces. The new class of tangent 
Levy models is called “discrete tangent Levy models”, as the jump sizes of the logarithm of the tangent process are 
restricted to finitely many values. As a result, for each fixed maturity, the corresponding Levy measure is purely 
atomic and can be represented by a finite number of parameters. There are several benefits of the new setting as 
opposed to the one considered in Section]^ On a theoretical level, the new drift restriction is simplified to a sum, 
as opposed to an integral which needs to be approximated numerically. In addition the mapping between option 
prices and Levy measure is simplified in this case. The latter makes the calibration problem somewhat easier, 
which, in turn, allows us to use non-parametric calibration procedure, which can potentially improve the quality 
of static fitting. The downside of this method is that, despite the existence of an explicit formula that connects the 
Levy density and options prices, solving the optimization problem associated with the non-parametric calibration 
is still very computationally expensive. We show how to deal with this problem further in the section. It is worth 
mentioning that discrete-space versions of tangent Levy models have been considered in and ll22ll . The results 
of lf3^ are limited to a single maturity, while the present setting includes multiple maturities. The theoretical results 
of II 22 II are very similar to the ones presented in this subsection. However, our choice of a convenient subclass of 
discrete tangent models and its subsequent numerical implementation are different. 

Similar to Section]^ herein, we assume that the true dynamics of S under the pricing measure Q are given by: 

St = So + f f Su-{e^ — l)[M{dx,du) — K{dx,du)], (3.1) 

Jo Jr 

where M is an integer-valued random measure with predictable compensator K{dx, du). The only difference be¬ 
tween •HD and (j^ is that, in the above expression, the compensator K{dx, du) may not be absolutely continuous 
with respect to dxdu. In fact, herein, we restrict our analysis to the compensators K of the form: 

N 

K{dx,du) = K^S^j {dx)du, (3.2) 


where every is a nonnegative predictable process, such that 

E f K^du < 00 , 

Jo 

^Please note that we cannot use 12.12) to calculate option prices, because, even though the calibrated Levy densities {rjt(T, •)} are double 
exponential, there is no reason to believe that the simulated r^’s remain double exponential. 
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and {x^, • • • , x^} is a finite subset of M which does not change with time. Clearly, x^’s correspond to the jump 
sizes and K^’s - to their intensities. At any fixed time t, a tangent model for the underlying is given by 

St = St + f ( Su-{e.^ — ^)[Nt{dx,du) — Kt{dx,du)] , (3.3) 

J t Jr 


for T G [0, T], where Nt{dx, du) is a Poisson random measure associated with the jumps of log S whose com¬ 
pensator is given by a deterministic measure Kt{du, dx). Of course, to be consistent with ( |3.2| l, we only consider 
Levy measures Kt of the form; 

N 

Kt{dx,du) = ''^^K{{u)6^j{dx)du, (3.4) 

f=i 

where every (•) is a continuous deterministic function. With a slight abuse of notation, we denote the collection 
of k{, for j = 1,..., N, by Kt- As before, we denote by (T, x) the call prices produced by the above tangent 

model at time t (cf. ( |2.3| l). The concept of a tangent model, then, requires that ( |2.4| l holds for all every t G [0, T). 
Next, we define the joint dynamics of St and Kt by 


St = So+ /o 4 Su- (e" - 1) Midx, du) - Ef=i idx)d 


4 (T) = k^T) + 4 ai{T)du + EZi fo j = l,...,N, 


(3.5) 


where B — (B^, ■ ■ ■ , B^) is an m-dimensional standard Brownian motion; (Kt)tG[o,T] ^nd (aDt^io^r] Pro¬ 
gressively measurable stochastic processes taking values in Bd (as defined in ( |6.7| l); and, for each j G {I, - , N} 

and n G {1, • • • , m}, (/3(’")tg[o f] ^ progressively measurable square integrable stochastic process with values 
in Hd (as defined in ( |6.8| l). 

As before, given the dynamics of {St, Kt), we need to ensure that they satisfy certain “consistency conditions” 
in order to produce models which are, indeed, tangent to the true model almost surely at all times. First, we make 
the following assumption. 

Assumption 2. The jump sizes {x^, ■ ■ ■ , are regularly spaced and have 0 at the center. In other words, the 
number of feasible jump sizes N is odd, and x^ = {j — M)dx for j = 1, - ■ ■ ,N, with M = {N l)/2 and let 
dx > 0 being the spacing between two feasible jump sizes. 

Note that x^ corresponds to jump size 0, therefore, the value of k^^ (T) can be set arbitrarily without changing 
the options’ prices. For convenience, we define it as 


k^{T) = -Y,kI{T), (3.6) 


Y.4{T)=0. (3.7) 

i=i 

Assumption 3. Almost surely, for every n = 1,..., m and almost every f > 0, we have: /3f= 0, Vfc < M/2 or 
k > 3M/2, with M ={N+ l)/2. 

Finally, we are ready to formulate the consistency conditions for discrete tangent Levy models. The theorem 
presented below, essentially, follows from Proposition 9 in ll22ll . up to some technical differences in the definition 
of the code-book dynamics ( |3.5| l. However, since the results of ll22l have not appeared in a publication, we present 
an alternative proof of this theorem in Appendix C. 

Theorem 2. Let Assumptions^and^hold and assume, in addition, that iSt)t^^Qf] martingale and 

Kt{T, x) > 0, almost surely for all t S [0, T) and almost all {T, x) G [t, T] x K. Then the processes {St, 
satisfying o and (O- are consistent, in the sense that ( |2.4| l holds, if and only if the following conditions hold 
almost surely for almost every t G [0, T) and all T G {t, T]: 
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7. Drift restriction: 


\{T) 


m .T 

^ ^ /3f’"(r) / j = i,...,7V. 

n=l w+l ^,,^ 3(iv+l) 

4 4 


(3.8) 


2. Compensator specification: 

3. Symmetry condition: 


Ki = kI (t), for all j 7 ^ M. 


X! - 1) = 0, n = l,...,m. 


(3.9) 


(3.10) 


w + l 3(JV + 1) 

4 —J — 4 


Notice that Theoremand equations ( |3.5| l provide a method for constructing a market-based model for call 
prices, or consistent model for implied volatility surface. We start by choosing {Pi ..., /3™)tg[o t]’ '^hat they 
are consistent with historical evolution, satisfy Assumption|^and the following linear constraints: 


E 


pr{T) = o, 


E 


Pl’^{T){e^ - 1 ) = 0 . 


-<j< 


3(iV + l) 


-<J< 


3(-W + l) 


Then, we use the drift restriction p.8[ ) and the second equation in p.5| l to generate the paths of (Kt)tg[o T]- Notice 
that p.7|i is satisfied with such choice of a and /3, which follows from 


N 


N m 


Y.al{T) = -J2J2 E 


\T) 


i=i 


j = l n=l N + 3(7V + lj 

4 4 


= -E E 




n=l iv+l^,.^3(JV+l) 

J 


where the last equality holds because the range of u covers regardless of the value of k. Having 

simulated the paths of t]’ compute the associated call prices Cl’^*. There is no need to simulate St, 

because, as discussed at the end of Subsection 2.1 for any value of St, the implied volatility in the log-moneyness 
variable can be computed by the inverse Black-Scholes formula applied to Cl’'^*. 


Similar to the discussion in Subsection 2.2 the implementation of discrete tangent Levy models consists of 
two stages - estimation and simulation. The estimation, in turn, consists of two parts: static fitting and dynamic 
fitting. All these steps are addressed in the following subsections. The market data consists of preprocessed SPX 
call option prices described in Subsection|2.3| 


3.2 Static fitting 

Recall that static fitting refers to the process by which we identify the additive measures Kt- This is done “statically” 
for each day t. As shown further in this subsection, in the case of discrete tangent Levy models, the PIDE for call 
prices can be simplified to a finite system of ordinary differential equations (ODEs). Moreover, the discrete tangent 
Levy models allow us to perform non-parametric calibration of tangent processes to market data, which helps us 
better understand the true structure of jump measures. By employing non-parametric htting with more variables, 
we allow more flexibility in the model hence we can obtain better fitting results. However, in order to obtain 
reasonable and consistent jump measures, a few technical difficulties associated with non-convex optimization 
problems and their potential instability need to be addressed. In this subsection, we explain how we overcome 
these difficulties in detail. 


We first describe the preliminary constructions. As specihed in Section 3.1 we assume the set of feasible jump 
sizes is equally spaced with zero at the center, so each can be expressed as (recall that M = is 


the center point with x^ = 0) 


—A + {j — l)dx with 


- 1 

A = {M — l)dx = — - — dx, 
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and {—A, A) is the truncation of the x-domain. The choice of N and dx requires some consideration. Ideally, we 
would like N to be as large as possible, to fully cover the structure of jumps, but a trade-off with the computation 
complexity and the limited data size is also important. As for the jump size spacing dx, it not only determines the 
value of A (the truncation of the x-domain), but also has a big impact on the jump intensities k^ ’s. On the one 
hand, dx should be large enough, so that {—A, A) covers the possible range of jump sizes and option prices. On 
the other hand, if dx is too large, the grid of jump sizes is too sparse, and the calibrated k might not reflect the true 
jump structure. A quick test on historical S&P 500 prices in the year of 2007 shows that the magnitude of most 
daily jumps falls between -0.035 and 0.029. Taking these facts into account, we decided to take N = 301 and 
dx = 0.005 after a few trials, so the center point M = 151 and the space truncation A — 0.75. However, the large 
number of variables makes the calibration (i.e. static fitting) nearly infeasible, as the procedure becomes extremely 
slow and the results are unstable. To address these issues, we reduce the number of parameters from N = 301 
to Nvar — 24, by assuming = 0 for large positive or negative jumps, and by dividing the jump sizes into 
small groups and making adjacent jump sizes share the same value of . The reduction speeds up the procedure 
significantly, and the calibration results look much more reasonable as the dimension of the optimization problem 
is greatly reduced. At the same time, the original grid with N = 301 and dx = 0.005 is kept sufficiently large, to 
allow for a good quality of fit to the market call prices. 

To perform the calibration, we need to develop a procedure to compute call prices from the calibration variables 
rI’s. Equation ( |8.3| l in Appendix C provides the following PDE for call prices: 

OrCf-"* (T, x) = ^ {T)e^^ (T, x - x^) - (T, x) • 4 (T) 

j:j^M 

+ [d, - l]C'f‘’"‘(T,x) • niT - t), (3.11) 

where H is the symmetry index defined by 

n{T - t) = - l) (3.12) 

k^M 


It is not hard to see that, if ( |3.8| l, ( |3.10| l and Assumption |^hold, then il( ) is a deterministic function. Notice that 
( |3.11| ) is a PDE because of the presence of the derivative {T, x) on the right-hand. Let us now perform a 

change of variables that eliminates this derivative and reduces the equation to a (multi-dimensional) ODE. Namely, 
we define 


Ct{T,x) :=C'r 


T, X— / n{u — t)du 


(3.13) 


This change of variables translates (|3.1 l|l to the following ODE of Ct{T, x): 


dTCt{T,x) 


^ 4iT)e-'CtiT,x-x4-CtiT,x) ^ 4iT) - ^iT - t) ■ Ct{T,x), 

j:j^M 


(3.14) 


which doesn’t contain any derivatives in the right hand side. Still, we will not use ( |3.14| i directly in computing 
option prices. Similar to Subsection 2.4 we would work with option time values instead of option prices. In our 
problem, we define the modified time value function as 


V^T) = CtiT, x^ + log St) - Stil - e"' )+. 


(3.15) 


The evolution of V at the grid-points is, then, given by 


'dTV:{T)= 4iT)e^'VtiT,x^-x4-VtiT,x^) ^ 4iT) 

- il{T - t) (Vt{T,x^) + St{l - e^')+) , i = ,N 
. V4t)=0, i = 


(3.16) 
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Let {Ti, T 2 , • • • ,Tl} denote the set of maturities available. We now describe the procedure to compute 

given Kt, using ( |3.16| l. Consider the time average of (T) between each pair of consecutive maturities Tj) 

(assuming Tq = t): 


m) 


1 


Ti - Ti_, 



{u)du, 


1<1 < L. 


(3.17) 


It is easy to see that the option time value doesn’t change if we substitute every Kj(-) to a piece-wise constant 
function, whose value between each T/_i and T; are given by Oj (Ti). With this substitution, the coefficients in the 
right hand side of the PDE ( |3.16| l become constant for T G (Ti_i, Tj). Hence, we can solve it to obtain an iterative 
formula for the solution at maturities T;: 


Vt{Ti) = _ i)G-\Ti)btiTi), Z = 1, • • • ,L, (3.18) 

where the matrix Gt(Ti) = has the following entries: 

^ B^Ti) + - f)] , 

and &t(Ti) = {bi{Ti),--- ,bN{Ti)Y is given by: 




.k-.k^M 


St eUTi)e^ 

k-.k^M 


^{1- 


with ^ 

n{Ti-t) = -—— f n{u-t)du. (3.19) 

j-l — JTi-i 

It is easy to see that Gt{Ti) and bt{Ti) only depend on {9l as does Vt{Ti). Notice that the true time value 
(the value by which the option price exceeds the intrinsic value) can be recovered from the modified time value by 






t)du -f log — S'* ^1 — e 


x^—J^ Q{u—t)du' 


(3.20) 


= Vt (T, x^) -f S'* (1 - n{u-t)du^ + ^ 

Thus, ( |3.18| l and ( |3.20| l establish the relation between the time values and the average jump intensi¬ 

ties 9{ (T;)’s. 

Let {Ti) denote the market time value at time t of the call option with strike Ste^^~^G-t)du 

maturity Ti. If any of these strike values are not traded in the market, we use linear interpolation in log-strike to 
obtain the value of (T/)’s. Increasing I, we formulate the calibration problem for each T/ separately: find 

nonnegative {9^ {Ti)}j^M which minimize the difference between the market time values {Vj™ {Ti)}^=i ™d 
the model time values under the symmetry condition 

^ - 1 ) =H(r,-f), (3.21) 

k^M 


with some fixed deterministic function ft. This symmetry condition ensures that ( |3.12| i holds, which, in turn, im¬ 
plies that p.8| l, p.l0[ l and Assumption[^hold. Notice that we formulate the calibration problem for each maturity 
separately because we can construct the time values {Ti)}S^^ using only {9{(Ti)}j^M and the already 

constructed (calibrated) time values for previous maturity, The preliminary optimization 

problem becomes: 


min 

Bt(Ti)>0 


N 

3=1 


(TO-Pr ’''(Tz)l' 


s.t. ^0^(TO(e"^-l) = H(r,-f), 

3^M 


(3.22) 
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where ujj = .^bid,^yask,j ,j. ^6 the weights associated with liquidities. As before, we put more weights 

It \i/ t \t/| 

on contracts with better liquidity, which is reflected by smaller bid-ask spread. Once the average intensities 9l (Ti) 
are constructed for all j and I, we recover |Kj (•) | using p.l7| l and interpolation across maturities. 

It is easy to see that ( |3.22| l is a non-convex optimization problem, and numerical algorithm may get stuck 
in a local minimum. In our implementation, we took three actions: (1) choosing an appropriate normalization 
function and working with normalized jump intensities, (2) carefully choosing the initial point, and (3) adding 
several penalization terms to regularize the problem. We now describe the three actions in detail: 

Action 1 Notice that the average jump intensities 9^ (T)’s can vary in size significantly across j (e.g. recall that 
the small jumps occur much more frequently than the large ones). This, in turn, may cause numerical 
difficulties in the optimization. To address this problem, we will use a normalization function denoted by 
Pi^) = P’ {x), and perform the optimization for the normalized average jump intensities 

9l{Ti) := 9l{Ti)/(P 

We choose the normalization function p to be given by 9q{Ti), which is obtained by solving p.22| l for the 
first maturity and the first day in the calibration period, without a symmetry constraint. This function is 
shown in Figure|^a. We then require the optimization variables (7/)}^ ^ to stay in a reasonable scale, to 
avoid the overflow problems. 

Action 2 Motivated by the choice of the normalization function, on each calibration day, we set the initial value 
for the optimization variable 9t{Ti) to be a vector of I’s. We use a calibrated 9t{Ti-i) as the initial value 
for 9tiTi). 

Action 3 To improve the convexity properties of the objective function and, therefore, make the optimizer be more 
likely to converge, we add the three penalization terms: 

Penalization term 1 


with 9{{Tq) = 1 
maturities T. 

Penalization term 2 

( M-l N \ 

Y,mT)-~9r\T)fQ{x^)+ Y. {0i{T)-9i+\T))^eix^)\ (3.24) 

j = l j=M+l J 


Fi{9t{Ti)) = Y (3.23) 

j^M 

. This penalization term also ensures that 9l{Tys is not varying too much across 


This penalization term also ensures that 9^ (T)’s is not varying too much across j. Here, 0(a:) is a fixed 
weight function, with which we penalize the differences at larger jumps more heavily - it is shown in 
Figure l^b. 

Penalization term 3 

^ 1 

This penalization term also keeps 9^5 away from zero, so that, if we simulate the future values of 0^, 
starting from u = t, they are less likely to touch zero. 

Summing up the above, we formulate the resulting optimization problem. For each day t G [0, T] and each 
maturity Ti, we run the penalized optimization 


min 

et(Ti)>o 


N 

ySt,pet,j 

i=i 



" + eiFl{9tiTi)) + €2F2{9t{Ti)) + e^F^Ti)) 


s.t. 


(3.26) 


16 












The coefficients of the penalization terms ei, £2 and £3 can be determined by a bisection method as follows. 


Algorithm 2: Algorithm for determining £ 1 , £2 and £3 

1 Run the unpenalized optimization ( |3.22[ ) without the symmetry constraint, and record the value of the 
objective function as /□; 

2 for 2 = 1 : 3 do 

3 £i = 5; 

4 / = lOO/o; 

5 while / > I.O 5/0 do 

6 Run the penalized optimization ( |3.26| l with the i-th penalization term only, and record the value of 

as /; 

Having run the Algorithm for different t and T;, we choose the following values: £1 = 0.3125, £2 = 
0.3125, and £3 = 0.0012. 


The average symmetry index function H( ) (introduced in ( |3.19| l) is then determined using t = 0 - the first day of 
the calibration horizon - by running the penalized optimization (3.26|l without the symmetry constraint and setting 


N 




(3.27) 


i=i 


We interpolate linearly to obtain the values of H(r) for all r G K+. Before we present the results, we summarize 
the static fitting procedure as follows. 

Algorithm 3: Algorithm for fitting 9 

1 Preprocess the market data according to ( |2.8| i; 

2 Determine £ 1 , £2 and £3 - the coefficients of the penalization terms - by Algorithm 

3 On t = 0, starting from the initial point dinit, specified in Action 2, run the penalized optimization ( |3.26| ), 
without the symmetry constraint, and compute the average symmetry index Cl via ( |3.27| ) and linear 
interpolation; 

4 For all days t G (0, T], run the constrained penalized optimization ( |3.26| l with 9init and Cl. Save the time 
series (0t)tG[o,T]- 

On Jan. 3, 2007, the first day of the calibration period, we fit 9q{Ti) to the prices of call options of six different 
maturities Ti, by solving ( |3.26| l without symmetry condition. Multiplying the results by p, we obtain the average 
jump intensities 0o(Ti)’s and the corresponding calibrated time values. With the family of 0o(Ti)’s, calibrated on 
the first day, we compute the average symmetry index Cl via \3.21) and linear interpolation, as shown in Figurej^c. 
Using the symmetry index, we run the constrained optimization ( |3.261 l for the subsequent days. The calibration 
results for the next day, Jan. 4, 2007, are shown in Figures [T^ - [^ It is clear that the fitting is very accurate 
in terms of matching the time values. It might be a little surprising to see the bimodal or even trimodal shape of 
the jump intensities; however, these were also observed in El. We can see that a discrete tangent Levy measure 
provides more flexibility and serves as a good nonparametric alternative to the parametric calibration performed in 
the previous section. 

Finally, we need to recover /tt’s from the calibrated 0(’s, via p.l7| i. To do this, we assume that each 

= Ki{T)/p> 

is exponential between two consecutive maturities T;_i and T;. In other words, we search for (T) in the follow¬ 
ing form: 

kiiT) = 

where clj, - ■ ■ , c^_i are constants. This assumption, along with (|3.17|i, leads to 


9llTi){Ti-Ti_,) 


4-1 



(T,-Ti_i) 



(3.28) 
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The above equation suggests that, as long as we have the value of kl (Tq) at one point Tg S [t, Ti], we will be able 
to back out k{ (T) for all T G \t, Ti] from {9{ (Ti)}[l^’s. Hence, for a fixed j, we will make a “guess” of kj (t) by 

1=1 

s.t. - 300 <4 <20, 1 = 1, L, (3.29) 

where 9l'^(Ti) is the T/]-average of the k^ (T), as given by ( |3.28| l, assuming that kj (t) = y. Once we obtain 

the optimal y*, we set kl (f) = y* and solve for c/’s via p.28[ ). The latter, in turn, gives us the values of kj (T) for 
all T G [t, T], Repeating the process for all j 7 ^ M and all t G [0, T], we obtain the time series of actual jump 
intensities An example of calibrated Kt is shown in Figure |T4^ a. 


running a single-variable optimization 


mm 

y>0 


3.3 Dynamic fitting 

Herein, we estimate the dynamics of Kt, given the historical sample obtained from option prices in the previous 
subsection. Recall that the dynamics of Kt are given by the second equation in ( |3.5| ). In the presence of the drift 
condition ( |3.8| l, the problem of dynamic fitting boils down to fixing the number of necessary factors m and finding 
the factors As in Subsection 2.5 we will use PCA to complete this task. Recall that PCA has to be 

applied to a stationary time series, hence we change variable from T to t = T — t and introduce 

Kt(r) = Kt{t + t) 


This gives us 

m 

dklir) = d^(r)dt + (3.30) 

n—1 

with 

(t) = a^ti* + t) + — + t). (3.31) 

For each day t, we choose 4 maturities, and for each maturity T, we have 24 jump intensities k^ (T)’s. So there 
are 4 x 24 = 96 points on the random surface Ak^ (T) = kl (T) — kl_i(T), which is quite large compared to the 
number of observations, 251. This could lead to inaccuracy in estimating the covariance matrix, thus hampering 
the estimation of /3"’s. To reduce the number of points, we discard the jump sizes whose jump intensities k^ are 
consistently very small. This gives us a much smaller set of 7 jump sizes, or 28 points on the surface of Akt- 
We perform PCA for the resulting time series {Awt}, to obtain the values of each /3^’”(T) at the chosen points 
{j, T). To obtain its values at other points, we simply interpolate /3 ’”(-) linearly. The results of PCA are shown in 
Figure [T 3 I We can see that the first 4 eigenvalues account for over 86 % of the total variance. Hence, we set 

to = 4 and /3" = \/A^/„, n = 1,2,3,4, 

where A„ is the n-th largest eigenvalue, and /„ is the corresponding eigenfunction. The so-obtained /3"’s are also 
shown in Figure[^as functions of x and r. After transforming /3 to /3 via ( |3.31| l, we apply the drift condition ( |3.8| l 
to compute the drift term a that guarantees the absence of dynamic arbitrage. A plot of a is shown in Figure [T3| 
Finally, we use (|3.31|) once more to obtain a from a. 


3.4 Monte Carlo simulation of implied volatility surfaces 


Using the results of static and dynamic fitting, herein, we simulate the future paths of implied volatility surfaces 
that are free of arbitrage and are consistent with the present and historical observations. Similar to Section 2.6 we 


incorporate a scaling factor 7 ^ in (|3.30[), to ensure the non-negativity of 


K^’s: 


dklir) =7i2a'’(T)df-f (3.32) 

n—1 
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where 



min 


l(T)Ae 


(3.33) 


with e = le ® and r = 0.5. Our strategy is now similar to the previous section. Namely, we use the explicit Euler 
scheme for (|3.32|l to simulate future k’s: 


«t+At('r) = + 7t XI (3-34) 

n—1 


with Af being one day. For each sample path, we convert simulated k to Kt{T) = kt{T — t) and compute the 
average intensity 6 defined in ( |3.17| l by numerical integration. We then apply the iterative formula ( |3.18| l to the 
average intensity 9 to compute the time values and, in turn, option prices {T,x) for different strikes e® 
and maturities T on any given day t. Finally, we invert the Black-Scholes formula to obtain implied volatilities. 
We perform the simulation starting from Dec. 28, 2007. Moving five days forward, the simulated /ct’s and the 
corresponding implied volatility surfaces are shown in FiguresandWe see that the simulated /ct’s maintain 
similar shape as the k on Dec. 28, 2007, and the implied volatility surfaces are consistent with what we normally 
observe in the market. Using this approach, we can generate as many samples as needed, and use them as scenarios 
of future implied volatility surfaces in various applications, such as pricing forward-starting options or solving 
problems of risk management. We illustrate this idea in Section]^ 


4 Empirical analysis of the performance of tangent Levy models 

In this section, we discuss the importance of consistency in modeling derivatives prices. As we know, an investment 
manager’s portfolio or a trader’s trading book often contains multiple financial derivatives written on the same 
underlying. As a simple example, an equity trader might hold a calendar spread and a butterfly spread at the same 
time. To properly manage the risk, one needs to understand the joint dynamics of these derivatives, for which a 
consistent modeling framework is crucial. Tangent Fevy models (as any market-based model) are built to achieve 
this goal precisely. This is due to the fact that not only present but also historical information contained in the time 
series of options’ prices is used in the estimation of model dynamics. Classical stochastic volatility models cannot 
capture the historical evolution of options’ prices, henc,e there is a reason to believe that market-based models 
would lead to better performance in portfolio management. To show that tangent Fevy models do indeed work 
better, here, we test the two tangent Fevy models implemented in Sections]^ and [fusing the following portfolio 
choice problem. The results are, then, compared against one of the most popular volatility models in the industry 
- the Stochastic Alpha Beta Rho (SABR) model. 

4.1 The variance-minimizing portfolio choice problem 

This example is a simplified Markowitz-type portfolio optimization problem. Consider a portfolio manager who 
needs to decide how he/she should balance a portfolio of SPX options so that its risk is minimized. Among the 
many definitions of portfolio risk, we adopt the one used in the classic Markowitz problem (for example, see 
Section 6.6 of ESll ') - namely, the standard deviation of the portfolio return over a given (future) time period. 
Notice that this is not a typical Markowitz portfolio problem, given we are not considering the trade-off between 
return and risk as a typical Markowitz problem would. As a matter of fact, we would assume that the portfolio 
manager lives in a risk-neutral world, so that the expected return is normalized. We admit that lacking excess 
return might make the example less exciting, but it helps us compare the model performance in an apples-to-apples 
fashion. With the normalized return, there is no need to worry about the impact of different market views portfolio 
managers might build into the investment decisions. Of course, without such a trade-off, there is a trivial solution 
to the portfolio choice problem - do not invest at all, reducing the risk to zero. To make the problem non-trivial, we 
require that the value of the portfolio at the time when it is constructed must be equal to a fixed number M. Such 
a restriction is relevant if the manager makes profits off the commission, proportional to the size of the investment 
portfolio he/she manages. For example, an option market maker might want to know the optimal inventory so that 
he/she can adjust the quoting strategy accordingly to reach the portfolio composition with minimal inventory risk. 
Or, a broker dealer might need to know her optimal position in options over the next several days to meet the risk 
and capital requirements. 
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We now formulate this problem mathematically. Let us assume that there are n options with the same maturity 
Tbut with different strikes Ki,Kn in the portfolio. Let Cu{Ki) be the time-u price of the itii-struck option, and 
let Wi be the quantity of this option in the portfolio, with a negative representing to short-selling. The weights uji 
have to be determined at the initial time d. The portfolio value at any future time t is simply Vt = X]r=i 
and the return over a u-day period is = Vd+ulVd- For simplicity, we assume that the risk-free rate and the 
dividend yield are both zero, so the expectation of is simply 1. For a given u G (0,T), to determine the 
portfolio weights, we need to solve the following convex optimization problem: 

min E(i?„ - 1)^ = ^ min E(14 - 
s.t. Vd = uFCd = M, 


where M £ K, M > 0 is the initial value of the portfolio. This is equivalent to 

min uF (4.1) 

tjGR" 

S.t. Vd=uj^Cd = M, 


where A„ = E[(C£;+„ — CdF is the covariance matrix of the time-d + u options’ prices. It is easy to see that the 
closed-form solution to the quadratic optimization (|4.1|l is 


MK^Cd 

CjkFCd 


(4.2) 


Thus, as it is well-known, the key to solving this optimization problem is to estimate the covariance matrix A„. To 
do this, we compute the sample covariance matrix using the time-d + u option prices simulated under each model. 
Note that, to obtain a fair comparison, the parameters of each model are only estimated using the options data prior 
to day d. Then, given N samples of the time-d -f u options’ prices. 


CU) = 




(j) 


j = 1, 




the sample covariance matrix is estimated as 

1 ^ 

A. = - Cd){C^^^ - Cdf. (4.3) 

^ i=i 

Different models generate different simulated paths of options’ prices, which then lead to different optimal weights. 
Naturally, how these optimal weights perform in the real world serves as an indicator of the model consistency. 
To be more specific, a consistent model should be able to generate portfolios with smaller standard deviation in 
the returns. To estimate the standard deviation of portfolio returns, we define the hgure of merit Q as the average 
realized deviation of the portfolio return in the testing period, i.e. 


Q = 




1 

Ntest 


Nt^st 

E -1)^ 


k=l 


(4.4) 


where Ntest is the number of trials, and is the actual portfolio return (given by market data) over a u-day period, 
with initial day dk and with the optimal weights F, obtained by ( |4.2| l on day dk- Different trials correspond to 
different initial days dk, i.e. 

1 " 

= (4.5) 

Recall that, by assumption, the mean of should always be 1. To make this assumption be consistent with the 
data, we choose a relatively small time horizon u. 
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4.2 Simulation algorithms 

As mentioned in the previous subsection, to find the optimal portfolio, we need to estimate the covariance matrix 
using simulated option prices. In this section, we describe the simulation algorithms for each model. 

• Double exponential tangent Levy model. For this experiment, we need to simulate both the underlying 
process S and the non-homogeneous Levy density k. For the double exponential tangent Levy model, in 
particular, we need to complete the following two steps to move one step ahead from ttot + At: 

- Step 1: Simulate the underlying process by 


St+At = St exp{- / (e"" - x)dxdt + ^ J^}. 

Here Kt{t, x) is the additive density for immediate maturity T = t, Nt is the number of jumps during 
the {t, t + Af] period, which has a Poisson distribution with parameter XAt, where A = /^ Kt{t, x)dx, 
and Jfc’s are the jump sizes having the disturbution x)dx. Notice that we approximate the jump 

component of log S with a compound Poisson process, which is reasonable given that the jump activity 
is finite in our setting. 

- Step 2; Simulate the Levy density Kt+At via ( |2.35| l. 

Simulating u days ahead requires repeating the two steps u times. We can then use the Fourier transform 
methods, as described in Subsection |2.6| to calculate time-u option prices, and estimate the covariance 
matrix to obtain optimal weights. 

• Discrete tangent Levy model. In this case, we need to complete two similar steps to move from ttot + At: 

- Step 1; Simulate the underlying process by 


( N Nt \ 

- l)Kl{t)At + ^Jk ■ 

j=i k=i j 

In the above, K((f)’s are the jump intensities for immediate maturity T = t, Nt is the number of 
jumps during the {t, t + Af] period, which has a Poisson distribution with parameter XAt, where 
A = (^)’ the jump sizes having the distribution ^ (dx). 

- Step 2; Simulate the intensity Kt+At via ( |3.34| i. 

Repeating the two steps u times, we simulate the Levy measure k and the underlying S' in u days. We then 
use the iterative formula ( |3.18| l to calculate time-u time values, from which we deduce the time-rt option 
prices and estimate the covariance matrix, to obtain the optimal portfolio weights. 

• SABR model. The simulation based on SABR model is slightly easier. SABR model, as proposed by Hagan 
et al. in ca, describes the dynamics of the forward price F and the volatility a as follows: 

dFt = atFfdBl 

dat = lyatdBf, (4.6) 

where F and a are correlated through dB^dB^ = pdt. ifTSl provides the following asymptotic formula for 
the time-f implied volatility under the SABR model: 


(Tt[K, T,Ft,at) 


_ot_ 

log2 Ft/K + log" Ft/K} ' 

ft [ (1-^)^ 1 2-3p^ 2 

i [ 24 {FtK)(^-d) 4 (FtA:)(i-4)/2 24 ^ 



iT-t)Y 


(4.7) 
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where K is the strike value, T is the maturity, Ft is the current spot level, and z and x{z) are defined as 


log§, 

at K 


(4.8) 


As for the parameters’ values, ifTSll suggests that /3 can be fixed in advance and ll32]| verifies empirically that 
this is a reasonable assumption. In our example, we will use two values of /3: /3 = 1 and /3 = 0.7. /3 = 1 is 
probably the most natural choice for equity market as it mimics a log-nomal model most closely, and /3 = 0.7 
is widely used on trading desks as it provides better results for risk management. The other parameters - the 
current volatility at, the volatility of volatility v and the correlation p - will be calibrated to market prices 
by minimizing the sum of squared differences between the market call prices and those produced by the 
model, calculated with ( |4.7| i. With the parameters calibrated on the initial day dk, the forward price and the 
volatility can be simulated as follows: 


Ft+At = 


at+At = 


(4.9) 


where Bl and Bt are independent. The time-u implied volatilities and option prices can then be computed 
via with the simulated spot Ft+u and volatility at+u- 


4.3 Results of empirical analysis 


In this section, we will go through the test procedure in detail and present the test results for the following models: 

• Double exponential tangent Levy model (DETL). 

• Discrete tangent Levy model (DTL). 

• SABR model with /3 = 1. 

• SABR model with /? = 0.7. 

Each model will be run in two periods: (I) Jan. 2007- Aug. 2008 and (II) Jan. 2011 - Dec. 2012. Eor each 
period, we use the first year’s data as a training sample, to estimate the parameters of the tangent Levy models, and 
we use the rest of the data as the testing sample, to compute the figure of merit Q defined in ( |4.4| i. The division 
between training and testing samples is shown in Table Please note that we cut off the first period at August 
2008 to reduce the impact of the financial crisis. The tests will be run on a portfolio of call options and underlying 
- referred to as a “(C + S) portfolio” - with three, four and five strikes. In each case, we pick every other strike 
starting from the strike closest to the underlying spot value (in other words, closest to at-the-money) at the moment 
when the portfolio is constructed. We pick these options because their market prices are most accurate. Assuming 
the set of available strikes is Ki < K 2 < ... < ATn and the spot S satisfies < S < Ki, Table [^illustrates 
the strikes used in each case. 

Eor all portfolios, we use a simulation horizon of it = 8 days, and, at the time dk when the portfolio is 
constructed, the options have maturity of T = dfc -|- it -|- 30, so that their time-to-maturity becomes 30 days when 
the given simulation period ends. We also assume the budget constraint M = 1. In addition to the figure of merit 
Q, we also check the average predicted deviation defined as 


P = 



Nt„st 


k=l 


(4.10) 


where is the set of optimal weights obtained via ( |4.2| l on day dk- The difference between Q and P is another 
measure of the accuracy of a model’s prediction. Besides the predicted and realized deviation, one may be inter¬ 
ested in how much the optimal portfolio weights fluctuate across the initial days dk- To measure this fluctuation, 
we define the average quantity oscillation index K: 


K = 




k=l 


1 


(fc+1)* 



(4.11) 
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where wf* is the quantity of the itTi-struck option in the optimal portfolio constructed on day dk- 


4.3.1 Period I 


For Period I, the parameters estimation for the two tangent Levy models is described in Sections 1^ and 1^ re¬ 
spectively. Following the simulation algorithms outlined in Subsection 4.2 for every initial day dk in the testing 


sample, we simulating 500 sample paths for the underlying and the option prices, using each model, and starting 
with the actual prices observed on day dk- In the simulation with tangent Levy models, we use the drift a and 
volatility /? estimated from the training sample (SABR model does not allow for any use of past prices). Using the 
simulated prices, for each model, we calculate the average predicted deviation P, according to ( |4.10| i, and estimate 
the optimal portfolio weights uj^* via •ED- Using these weights, we construct the corresponding and follow its 
value using the actual market prices in the time period {dk^ dk -f m], to compute the average realized deviation Q 
via ( |4.4| i. 

The results are shown in Table It is easy to see that, for a portfolio with 5 strikes, the tangent Levy models 
produce much smaller values of Q than the one produced by a SABR model, indicating that the tangent Levy 
models do a much better job at finding the minimal-variance portfolio. This can also be seen in Figure [T^ which 
shows that the distribution of realized returns is much more concentrated around 1 under the tangent Levy models 
than under the SABR model. Furthermore, if we look at the difference between Q and P, we can see that it is 
much smaller for tangent Levy models than for the SABR model. This suggests that the tangent models produce 
a more reliable prediction of the risk of an option portfolio (as measured by the standard deviation of its return) 
than the SABR model. Besides a small return fluctuation, another nice feature of tangent Levy models is the 
stability of optimal option quantities across the initial days dk- Figure [^shows the optimal quantities of options 
and underlying index, in the portfolio with 5 strikes, across all initial days in the testing period, for every model. 
Similarly, Table|^shows the average quantity oscillation K (defined in ( |4.1 l| i) for all portfolios and all models. It 
is easy to see that the portfolio weights constructed via the tangent Levy models are much more stable than those 
constructed using the SABR model. This can be explained by the fact that the parameters of tangent Levy models 
are estimated from both the present and historical option prices, while the classical stochastic volatility models, 
such as SABR, can only be calibrated to the option prices available on day dk- It is well known (and obvious 
intuitively) that an estimate based on a larger sample is more robust. Thus, the ability of tangent Levy models to 
be fitted to the historical options prices makes their output (in this case, the optimal portfolio weights) more stable. 

Tables]^ and 1^ also show that the difference between the performance of tangent Levy models and the perfor¬ 
mance of SABR model shrinks as the number of strikes in the portfolio decreases. This is not a surprise: as the 
number of strikes decreases, the number of degrees of freedom in the dynamics of option prices, which have to 
be captured by the model, decreases as well. Eventually, for a very small number of strikes, the SABR models 
are doing relatively well. However, even in the case of 3 strikes, the tangent Levy models do at least as good as 
SABR (although at a higher computational cost). Thus, the real benefit of using tangent Levy models is, of course, 
only visible when the number of options in the portfolio is relatively large. Figure [^provides a visual explanation 
for tangent Levy models’ outperformance. It shows the simulated call option prices, as functions of strike, at the 
end of the simulation period in 500 sample paths under the double exponential tangent Levy model and under the 
SABR model with /3 = 1. It is easy to see that the SABR model allows for very limited shapes of call price curves, 
while the tangent Levy model is able to generate a lot more various shapes. It is the lack of variety of different 
scenarios for the joint evolution of call prices (not simply the lack of parameters in the model) that prohibits the 
classical stochastic volatility models, such as the SABR model, from capturing the true dynamics of option prices 
(or implied volatility surface) contained in the historical data. 

It is worth noting that the discrete tangent Levy model doesn’t perform as well as the double exponential tan¬ 
gent Levy model, even though it is still superior to SABR model (assuming sufficiently many strikes). This can be 
explained by the fact that the results of the nonparametric static fitting of a tangent Levy density, discussed in Sub- 
section [T^ often lack stability, introducing additional noise into the model and, hence, damaging its performance. 
It could be an interesting area for future work to improve the robustness of static fitting of discrete tangent Levy 
models. 


4.3.2 Period II 

Herein, we repeat the same analysis for Period 11. The main purpose of this analysis is to show that the outperfor¬ 
mance of tangent Levy models is not due to our choice of a testing period, but that it is a persistent property. First, 


23 







we need to estimate the parameters of the two tangent Levy models using the data of year 2011. The calibration 
procedures are exactly the same as the ones described in Sections|^and[^ so we only present the main results here. 

• In the case of double exponential tangent Levy model, the first three eigenmodes explain over 93% of the 
variance. The eigenvalues and the eigenmodes are shown in Figure [T^ and the corresponding drift term a 
is shown in Figure]^ a. Comparing to Figures]^ andwe see that these results are almost the same as for 
the year 2007, suggesting that this model is very robust. 

• In the case of discrete tangent Levy model, the first three eigenmodes explain over 86% of the variance. The 
eigenvalues, the eigenmodes and the drift term are shown in Figures [2T|and[20|b. Note that these graphs are 
slightly different from those we obtained for year 2007 (Figures[^and ??), although having similar patterns. 
The lack of robustness can be explained by the fact, in the case of DTL models, we use a non-parametric 
fit to the data, which may offer a better fit quality, but, at the same time, is known to be less robust than the 
parametric estimation. 


Once the estimation is completed, we can repeat the same simulation and testing procedures as in Subsection 
|4.3.1| to obtain the results shown in Tables and |7j as well as in Figures andThese results confirm the 
finding of Subsection 4.3.1 for sufficiently many strikes in the portfolio, the tangent Levy models do a much 


better job at finding a portfolio with smallest variance, their predictions for the variance are more reliable, and the 
portfolio weights are more stable. 


5 Conclusion 


In this paper, we implement and test two types of market-based models for European-type options. These types of 
models can be viewed as numerically tractable specifications of the tangent Levy models proposed in lU and lO. 
In particular, they provide a method for generating Monte Carlo samples of future implied volatility surfaces, in 
a way that is consistent with their past and present values. We estimate the parameters of these models using real 
market data, for two periods: 2007-2008 and 2011-2012. The estimation procedure is described in a lot of detail, 
so that it can be reproduced by any interested reader. 

In addition, we use the estimated models and the real market data to conduct an empirical study using, whose 
goal is to compare the performance of market-based models with the performance of classical stochastic volatility 
models. We choose the problem of minimal-variance portfolio choice as a measure of model performance and 
compare the two tangent Levy models to S ABR model. Our study demonstrates that the tangent Levy models do a 
much better job at finding a portfolio with smallest variance, their predictions for the variance are more reliable, and 
the portfolio weights are more stable. To the best of our knowledge, this is the first example of empirical analysis 
that provides a convincing evidence of the outperformance of the market-based models for European options using 
real market data. 

Our work is also subject to certain limitations, which suggest directions of future research. One of the most 
serious limitations of our work is the lack of numerical stability at the stage of static fitting, discussed at the end 


of Subsection 4.3.1 To address this issue, one has to come up with a family of tangent Levy densities that is rich 
enough - so that it can approximate well the option prices observed in the market - and, at the same time, not too 
large - so that the calibration procedure is numerically tractable and more stable. This is a balance that seems to be 
hard to find. Another extension is to search for other families of tangent models - not necessarily based on Levy 
processes. Eor example, it is not hard to find a combination of arbitrage-free prices of two call options, with the 
same maturity and different strikes, which cannot be approximated with an arbitrary precision (i.e. both prices at 
the same time) by the exponential Levy models. This means that, in principle, the market prices of call options 
(even with the bid-ask spread) may be such that there is no tangent Levy model that can match them (provided 
they contain more than one strikes). This, in turn, motivates the search for other families of models, which can 
always fit an arbitrary family of arbitrage-free option prices. An example of such family is provided in ||7], but the 
existence and description of consistent dynamics within this family of models remains an open question. 


6 Appendix A 

Here, we define the Banach spaces associated with tangent Levy processes. 


24 






• Bo is a Banach space of Borel measurable functions satisfying 


II/IIbo := / {\x\Al)\x\{l + e^)\f{x)\dx <oo. (6.1) 

JR 

• B is 3. Banach space of absolutely continuous functions / : [0, f ] — Bq satisfying 

II/IIb := II/(0)|1bo \\-^f{u)\\Bodu<oo. (6.2) 

• Ho is a Hilbert space of Borel measurable functions / : K —>■ K satisfying 

ll/ll«o:= / \x\\l + en^\fix)\^dx<oo. (6.3) 

JR 

• "H is a Hilbert space of absolutely continuous functions / : [0, T] —>■ I-Lq satisfying 

l!/ll« := ll/(0)ll«o \\-^f(.u)\\nJu<oo. (6.4) 

• (^([O, T]) is a Banach space of continuous functions / : [0, f] —>■ M satisfying 

sup |/(a;)| < oo. (6.5) 

£CGR 


• T]) is a Hilbert space of absolutely continuous functions / : [0, T] —)■ K satisfying 

l/(0)P+/ lJ^/(w)P<oo. (6.6) 

Jo du 

• Bd i^ ^ Banach space of absolutely continuous functions / : [0, T] —M satisfying 

II/IIb^ := |/(0)| \-^f{u)\du<oo. (6.7) 

Here the subscript d is used to indicate the “discrete” models. 

• Hd is the Hilbert space of absolutely continuous functions / : [0, T] —)■ K satisfying 

II/IIL ■= l/(0)P \-^f{u)\^du<co. (6.8) 

We know that Ho C Bq, H C B, C C'([0,T']) and Hd C Bd- In addition, it is not hard to see 

that the completion of Ho is Bq with respect to the norm || • Hgg. Similarly, the completion of H is B with respect 
to II • lie, the completion of T]) is C([0, T]) with respect to the “sup” norm, and the completion of Hd is 

Bd with respect to the II • IIb^ norm. Hence, we conclude that the couples ("H, B), (kk^’^([0, T]), C([0, T])), and 
{Hd, Bd) are all conditional Banach spaces (see III 5.3 in 1^ for definition). 


7 Appendix B 

Proof of Lemmaj^ The proof is similar to the one given in ifT^ except that Z{T) = fiT + follows a 

gamma distribution in the absence of the diffusion term. The tail probability is given by 


F{Z{T) > a} = T'(/r, A,p, Ai, A 2 ;o,r), 


25 


with 'I' given in ( |2.13| ). If we set Vi = exp(li) for z = 1, • • • , N, the drift term has to satisfy /r = — AE[yi — 1] for 
St to be a martingale, so the dynamics become 


dSt = 


-AE[yi - l]St-dt + St- 


■ d 


■ Nt 


. 1=1 


Let C = E[yi — 1] = ~ Using results on equivalence of measures for compound Poisson 

processes (see Proposition 9.6 in iflOll for example), we can see the time-f price of a call option with maturity T 
and strike K is 


CtiT,K) = E[(5t - K) + \Tt] = E[^Tls,>tf|J't] - E[Kls^>K\J^t], 

= Sf's(^-XC,X\p*,XlXhlog{^),T-t 

- ^ ^ ^-AC, X,p, Ai, A 2 ; log(^), T -t] , 

where A* = A(C + 1), P* = • xEt ’ •^1 = ^1 ~ A 2 = A 2 + 1. 


8 Appendix C 

Proof of Theorem]^ Let us introduce the double exponential tail function ip, for any Levy measure v (defined on 
M \ {0}), with finite exponential moments, and all x S K \ {0}: 


ip{v,x) := 


- e-^)v{dz) x<0 

J^(e^ — e^)iy(dz) x > 0. 


( 8 . 1 ) 


Notice that, since is a Levy measure, ip{v, •) S L^(]R). In addition, ip(y, ■) determines v uniquely. Using the 
above definition of ip, we repeat the derivations in 0 , to derive the associated partial integro-differential equation 
(PIDE) for call prices (T, x) in the (T, x) variables: 


5TCf*>"*(r,x) = J^iP{Kt{T),x- y) (9L - d.) Cf^^^\T,y)dy 


( 8 . 2 ) 


where we treat Kt{T) as a measure: Kt{T) = ■ Integrating by parts in ( |8.2| i, we make use of the 

discrete structure of Kt{T) to replace the integral with a summation and obtain the following PDE: 


+ [9,-l]Cf‘’''‘(T,x)- E Ki{T){e-^-l) (8.3) 

j-.jitM 

Cf^’^^{t,x)= {St-e-)+. 


Similar to a, we introduce Aj (T, x) = —dxCt’'^* (T, x). Differentiating and taking Eourier transform in : 
on both sides of (|8.3|l, we obtain 


5TAf‘’"‘(T,0= Af 


Af‘’''‘(f,0 = 

whose solution is given by 


E Kt(^) -l + (27rzC-l)(e"^"-1)) ) 


g(l- 27 ri?) log St 

l-27ri^ ’ 


^(1 —27ri^) log 5t ( . y \ \ 

Af*’''‘(r,^) = ——exp ( ^ Kl{u)du ■ - (1 - 27ri^)e"^' - 27ri^j j . (8.4) 
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Repeating the first part of the proof of Theorem 4.7 in ||4], we see that the dynamics of {St, Kt) are consistent with 
the true underlying dynamics if and only if 


\ 7b /tg[0,T) 

is a local martingale for any (j) G S, where S is the Schwartz space of fast-decaying functions (cf. ID). This 
condition, in turn, is equivalent to the drift term in the semimartingale decomposition of (T, being zero. 

An application of the generalized Ito’s formula (cf. Theorem III.5.4 in IItTI I to gives the following 

expression for the drift: 


Tt{T,0 :=Af‘’«‘(r,0 


^ Ki{t) - (1 - 27rzOe"' - 27rze) 


r-i^M 
rT 


+ ^ a{{u)du ■ 27rz5)a:^ — (1 — — 27 rf^^ 




^ m f \ 

+ T E E / Pr(.u)du ■ - (1 - - 27r*c) 

^u=i\r.j^Mdt 

+ X! / - e'^(l - 27ri^) - 27rz^^ . 


jttM 

The normalized drift term T t {T, ^ )/Af* {T, is absolutely continuous as a function of T G [t,T], therefore the 
condition ri(r, = 0 is equivalent to 


and 


A direct calculation shows that 


limTA rt(r,^) ^ ^ 

\imTitAf^’^^{T,0~ ’ 


(8.5) 


( 8 . 6 ) 


— e'^(l — 27rf^) — 27rz^)Ar/(52,j(dx) = —27r(27r^^ + f e j ;a; ] dx, 

7b \j^m ) 

K{(t) ~ (1 ~ 2Tri^)e^^ — 27rf^'j = —27r(27r^^ + *0 [ ( y ^ i'^ii^)dxG,x \ dx. 

dR J 

Along with the uniqueness of Fourier transform and its inverse, the above equation shows that ( |8.5[ l is equivalent 
to the compensator condition p.9| l. Next, ( |8.6| l gives us: 

aiiT) ■ - (1 - 27rf^)e“' - 27^^e) 


r-j^M 


= - E E - (1 - 27r*e)e" - 2ni0 


n —1 \k\k^M 


- (1 - 2mCje^' - 2m£,) 


(8.7) 


with /f(’"(T) = f](’^{u)du. The equation holds for all ^ S K. Dividing both sides by and letting ^—>- 00 , 

we obtain 

0 = E ( E - 1) I ( E - 1) 

n=l \k-.kt^M ) \ld=jtM 
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= E E - 1) VO < < < T < r. 

n=l \k-.k^M J 

from which we conclude that /3 has to satisfy the symmetry condition ( |3.10| ). This symmetry condition also allows 
us to simplify (|8.7|) to 


aiiT) ■ - (1 - 2Tri^)e^' - 2 ^^^) 

j:j^M 


= -E 


k-.k^M 



[ ^ -e^') 


Now, dividing both sides of the above equation by and let ^ > oo, we see that a also has to satisfy the symmetry 
condition 

N 

Y oPATW' - 1) = E - 1) = 0. (8.8) 

j:j^M k=l 

Recall that we defined so that p.7| l holds: J2f=i (^) = This relation is preserved for all t G [0, T] if and 
only if 

N N 

Y<^t{T)=0, ^^r(T) = 0, (8.9) 

i=i i=i 

Substituting p.8| l, p.l0| i and p.9| l into p.7| i, we obtain 

N m / N \ / ^ 

Y^t{T) ■ " “E 

3 = 1 n=l \k=l / \i=l 



Both right and left hand sides of the above equation can be expressed as linear combinations of }i=i, - ,n- 

The latter functions are linearly independent, therefore, p.l0| l is equivalent to a system of equations, in which we 
equate the coefficients in front of every basis function. This, in combination with Assumption yields p.8| ). 
Finally, we notice that p.l0|i and p.8[) imply p.8|l. As p.9|l always holds (by the definition of k^), we can reverse 


the above derivations to show that ( 3.10| l and p.8| l imply p.7| i. Thus, 
P.10[) and p.8|l, which completes the proof of the theorem. 


have shown that p.6[) is equivalent to 


9 Appendix D 


Table 1: Time periods 



Jan. 2007 - Aug. 2008 

Jan. 2011 - Dec. 2012 

# of days 

419 

502 

Range of SPX spot price 

$1214.9-$1565.2 

$1099.2-$1465.8 
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Figure 1: Symmetry index ^ as a function of time to maturity with double exponential model 


Table 2: Testing periods 


Period 

Training period 

Testing period 

I 

II 

Jan. 2007 - Dec. 2007 
Jan. 2011 - Dec. 2011 

Jan. 2008 - Aug. 2008 
Jan. 2012 - Dec. 2012 


Table 3: Strikes used in each portfolio 


# of strikes 

Strikes used 

5 

if, 3 (call), X,i(call), X,+i(call), iT,+ 3 (call), iT,+ 5 (call) 

4 

if, 3 (call), X,i(call), X,+i(call), iT,+ 3 (call) 

3 

i^,_i(call), iT.+i(call), iT.+ 3 (call) 


Table 4; Average deviation of (C + S) portfolio in period I 



# of strikes 

DETL 

DTE 

SABR (/3 = 1) 

SABR (/3 = 0.7) 

Averaged realized 
deviation Q 

5 

0.55% 

1.13% 

84.97% 

111.42% 

4 

0.54% 

1.07% 

4.69% 

24.43% 

3 

0.64% 

0.91% 

2.18% 

10.50% 

Averaged predicted 
deviation P 

5 

0.87% 

0.88% 

0.19% 

9.33% 

4 

0.88% 

0.92% 

0.30% 

9.66% 

3 

1.05% 

1.03% 

0.53% 

10.29% 


Table 5: Average quantity oscillation K (as defined in (411 1 ) in (C + S) portfolio in Period 1 


# of strikes 

DETL 

DTE 

SABR (/3 = 1) 

SABR (P = 0.7) 

5 

0.0039 

0.0082 

1.1747 

2.3846 

4 

0.0038 

0.0081 

0.1339 

0.4629 

3 

0.0027 

0.0075 

0.0263 

0.0807 


Table 6: Average deviation of (C + S) portfolio in Period II 



# of strikes 

DETL 

DTE 

SABR (P = 1) 

SABR (P = 0.7) 

Average realized 
deviation Q 

5 

0.41% 

0.63% 

9.07% 

33.22% 

4 

0.42% 

0.57% 

3.51% 

17.61% 

3 

0.42% 

0.60% 

0.90% 

5.22% 

Average predicted 
deviation P 

5 

0.79% 

0.48% 

0.36% 

7.98% 

4 

0.79% 

0.55% 

0.43% 

8.11% 

3 

0.94% 

0.65% 

0.62% 

8.46% 
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Table 7: Average quantity oscillation K (as defined in (|4. 1 l|i) of (C + S) portfolio with 5 strikes in Period II 


# of strikes 

DETL 

DTE 

SABR (/3 = 1) 

SABR (P = 0.7) 

5 

0.0011 

0.0020 

0.1410 

0.6642 

4 

0.0012 

0.0023 

0.0537 

0.2736 

3 

0.0011 

0.0021 

0.0145 

0.0474 


Levy Density r| (on 20070104 for maturity 20070120 ) 


Levy Density T| {on 20070104 for maturity 20070217 ) 




(a) Maturity Jan.20,2007 

Levy Density T| (on 20070104 for maturity 20070317 ) 


(b) Maturity Feb. 17,2007 

Levy Density r| (on 20070104 for maturity 20070616 ) 




(c) Maturity Mar. 17,2007 

Levy Density r| (on 20070104 for maturity 20070922 ) 


(d) Maturity Jun. 16,2007 

Levy Density r| (on 20070104 for maturity 20071222 ) 




(e) Maturity Sep.22,2007 


(f) Maturity Dec.22,2007 


Figure 2: Calibrated densities 77 for DETL model on the second day, Jan. 4, 2007 
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Time Value {$) Time Value {$) Time Value {$) 


Time Value (on 20070104 for maturity 20070120 ) 




Log-moneyness 


(a) Maturity Jan.20,2007 


(b) Maturity Feb. 17,2007 



Log-moneyness 


Time Value (on 20070104 for maturity 20070616 ) 



(c) Maturity Mar. 17,2007 (d) Maturity Jun. 16,2007 

Time Value (on 20070104 for maturity 20070922 ) Time Value (on 20070104 for maturity 20071222 ) 




(e) Maturity Sep.22,2007 (f) Maturity Dec.22,2007 

Figure 3: Calibrated time values for DETL model on the second day, Jan. 4, 2007 
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(a) T} and k, for small jumps 


(b) 7 } and k for large jumps 



600^ 



(c) Levy density 77 (d) Additive density k reconstructed from 77 in (c) 


Figure 4; Calculating Levy density k from rj 
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10000V, 



(a) Percentage of variance explained by the eigenmodes 


(b) The first eigenmode scaled by \/\\ 



(c) The second eigenmode scaled by \/\2 


(d) The third eigenmode scaled by \/\z 


Figure 5: Eigenvalues and eigenmodes of Aii for DETL model 



Eigure 6; The drift term a for DETL model 
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1200 . 0 . 24 . 



(a) Calibrated k. 


(b) Calibrated Implied Volatility Surface 


(c) Simulated (1st day) 


(d) Simulated Implied Volatility Surface (1st day) 





(e) Simulated k (2nd day) 


(f) Simulated Implied Volatility Surface (2nd day) 


Figure 7: Simulated k’s and implied volatility surfaces using DETL model (1) 
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(a) Simulated k (3rd day) 


(b) Simulated Implied Volatility Surface (3rd day) 



(c) Simulated tz (4th day) 


(d) Simulated Implied Volatility Surface (4th day) 




(e) Simulated k (5th day) 


(f) Simulated Implied Volatility Surface (5th day) 


Figure 8: Simulated k’s and implied volatility surfaces using DETL model (2) 
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Time Value (on 20070104 for maturity 20070120 ) 




(a) Maturity Jan.20,2007 


(b) Maturity Feb. 17,2007 




(c) Maturity Mar. 17,2007 


(d) Maturity Jun. 16,2007 



Time Value (on 20070104 for maturity 20071222 ) 



(e) Maturity Sep.22,2007 


(f) Maturity Dec.22,2007 


Figure 11; Calibrated time values for DTL model on the second day, Jan. 4, 2007 
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(a) The first eigenmode scaled by v^Ai 


(b) The second eigenmode scaled by \/\2 



Time To Maturity (c) 



Time To Maturity (t) 


(c) The third eigenmode scaled by a/As (d) The fourth eigenmode scaled by v^Aj 

Figure 12: Eigenvalues and eigenmodes of A/i for DTE model 



X to' 



(a) Percentage of variance explained by the eigenmodes (b) The drift term ct for DTL model 

Figure 13: Eigenvalues and eigenmodes of Ak for DTL model 
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70- 

60„ 


50- 





(c) Simulated (1st day) 


(d) Simulated Implied Volatility Surface (1st day) 




(e) Simulated k (2nd day) 


(f) Simulated Implied Volatility Surface (2nd day) 


Figure 14: Simulated k’s and implied volatility surfaces using DTL model (1) 
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250 . 



0 . 35 . 


I" 0.3. 
I 0.25. 



(a) Simulated k (3rd day) 


(b) Simulated Implied Volatility Surface (3rd day) 


250. 



0.35. 


S' 0.3. 




(e) Simulated k (5th day) 


(f) Simulated Implied Volatility Surface (5th day) 


Figure 15: Simulated k’s and implied volatility surfaces using DTL model (2) 
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Return distribution of the calls + underlying portfolio under double exponential tangent Levy 


90- 

80- 

70- 

60 - 
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§ 50 - 
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40- 

30 - 

20 - 

to¬ 
ol-1--1-1- B-1-1-1-1- 

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 

Return 


Return distribution of the calls -i- underlying portfolio under discrete tangent Levy 


0.6 0.8 


(a) Under double exponential tangent Levy model 


(b) Under discrete tangent Levy model 




(c) Under SABR model with ^ — 1 (d) Under SABR model with ^ = 0.7 

Figure 16: Distribution of the 8-day returns of (C -i- S) portfolio with 5 strikes in Period I. Different scales are used 
to show more details. 
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(a) Under double exponential tangent Levy model 


(b) Under discrete tangent Levy model 




(c) Under SABR model with ^ — 1 (d) Under SABR model with /3 = 0.7 

Figure 17; Option quantities in (C + S) portfolio with 5 strikes in Period I. Different scales are used to show more 
details 



(a) Under double exponential tangent Levy model 


(b) Under SABR model with (3—1 


Figure 18; Terminal option prices in (C + S) portfolio, as functions of strike, simulated using 500 sample paths 
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(a) Percentage of variance explained by the first eigenmodes 


(c) The second eigenmode scaled by -\/A^ 


(d) The third eigenmode scaled by y /^3 


Figure 19: Eigenvalues and eigenmodes of Ak under DETL, estimated using 2011 data 






(a) Under double exponential tangent Levy model 


(b) Under discrete tangent Levy model 


Eigure 20: The drift terms a under DETL and DTE, estimated using 2011 data 
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Factor p, (i,x) 




Time To Maturity (t) 


(a) Percentage of variance explained by the first eigenmodes 


(b) The first eigenmode scaled by \/\\ 



(c) The second eigenmode scaled by -\/A^ (d) The third eigenmode scaled by \/\z 

Figure 21: Eigenvalues and eigenmodes under DTE, estimated using 2011 data 
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Return distribution of the calls + underlying portfolio under discrete tangent Levy 


(a) Under double exponential tangent Levy model 


(b) Under discrete tangent Levy model 



.6 1.S 2 



(c) Under SABR model with 0—1 (d) Under SABR model with 0 — 0.7 

Figure 22: Distribution of the 8-day returns of (C -i- S) portfolio with 5 strikes in Period II 
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(a) Under double exponential tangent Levy model 


(b) Under discrete tangent Levy model 




(c) Under SABR model with ^ — 1 (d) Under SABR model with ^ = 0.7 

Figure 23: Option quantities in (C + S) portfolio with 5 strikes in Period 11. Different scales are used to show more 
details 
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